cond-mat/0405363 



in 
o 
o 

(N 
X> 

PL, 
00 



Propagation of Ripples in Monte Carlo Models 
of Sputter Induced Surface Morphology 
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Periodic ripples generated from the off-normal incidence ion beam bombardment of solid sur- 
faces have been observed to propagate with a dispersion in the velocity. We investigate this ripple 
behaviour by means of a Monte Carlo model of the erosion process, in conjuction with one of 
two different surface diffusion mechanisms, representative of two different classes of materials; one 
is a Arrhenius-type Monte Carlo method including a term (possibly zero) that accounts for the 
Schwoebel effect, the other a thermodynamic mechanism without the Schwoebel effect. We find 
that the behavior of the ripple velocity and wavelength depends on the sputtering timescale, qual- 
itatively consistent with experiments. Futhermore, we observe a strong temperature dependance 
of the ripple velocity, calling for experiments at different temperatures. Also, we observe that the 
ripple velocity vanishes ahead of the periodic ripple pattern. 

PACS numbers: 05.10.-a,68.35.-p,79.20.-m 
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I. INTRODUCTION 

There has been much scientific activity for quite some 
time now, on the features of surface morphology result- 
ing from the bombardment of a solid surface by a col- 
limated beam of intermediate energy ions, at normal 
and oblique incidence to the solid surface 0, Q- The 
phenomenon is an essential constituent of several sur- 
face analysis, processing and fabrication techniques, such 
as ion beam aided deposition, surface catalysis, sputter 
cleaning, etching and deposition. 

Normal incidence ion bombardment of non-metallic 
substrates often results in an interlocking grid of hillocks 
and depressions, which have been demonstrated to be an 
attractive alternative to the spontaneous growth of self- 
organized quantum dots on semiconductor surfaces in the 
Stranski-Krastanov growth mode Off- normal inci- 
dence ion bombardment of such non- metallic substrates, 
however, gives rise to the formation of quasiperiodic rip- 
ples [IIIE 01111113,113 with orientation that depends 
on the angle of incidence of the ion beam. For incidence 
angles less than a critical angle, 6 C |l2j, the wavevector 
of the ripples is parallel to the projection of the ion beam 
direction on the surface plane while for incidence angles 
greater than 8 C , the wavevector of these ripples is oriented 
perpendicular to the projection of the ion beam direc- 
tion on the surface plane. On the other hand, ripples are 
observed on metallic substrates at normal incidence ion 
bombardment, and these ripples are rotated by changing 
the subtrate temperature [lj, 0, ; a probable conse- 
quence of the symmetry-breaking anisotropy in surface 
diffusion. The wavelengths of the observed ripples, in all 
cases, is of the order of tenths of micrometers. 

However, a number of experimental studies [TtilrillTsl 
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IX9| have demonstrated that under certain ion bombard- 
ment conditions, ripples are not formed; the surface un- 
dergoes kinetic roughening with interesting scaling prop- 
erties. All these observations point to the possibility of 
several phases in the surface topography evolution, with 
phase boundaries defined by the bombardment condi- 
tions, and with little or no dependence on the material 
composition, surface chemistry, defects or chemical re- 
actions on the surface. These features are understood, 
from insightful theoretical descriptions 0, I20I l2l) , as be- 
ing governed by the interplay and competition between 
the dynamics of surface roughening on the one hand and 
material transport during surface migration on the other. 
Ion bombardment tends to roughen the surface, while 
surface diffusion leads, in general, to surface relaxation 
For sufficiently low ion energies, the sputtering phe- 
nomenon is the dominating mechanism |2J. However, if 
the flux is low at such energies, then the enhanced de- 
fect mobility can result to domination by surface diffu- 
sion which may cause the overall scaling behavior of the 
surface profile to be uniquely determined by the nonequi- 
librium biased diffusion current, independently of the mi- 
croscopic origin p^ . 

Recently surface ripples generated during Gallium ion 
beam erosion of Silicon were observed to propagate with 
a ripple velocity that scales with the ripple wavelength 
as v ~ X k , where k ~ initially, and k — —1.5 after a 
crossover wavelength A c ~ lOOnm [23(. This velocity dis- 
persion has been ascribed to an indication of a continuous 
transition to a rising non-linear contribution in surface 
erosion H [2j|. Motivated by this experimental result, 
we study ripple propagation by means of a recently in- 
troduced, discrete (2+1) dimensional Monte Carlo (MC) 
model 24] of the sputtering process, and two different 
solid on solid models of surface diffusion; for details see 
below . We focus on intermediate times, where the transi- 
tion from linear to non-linear regimes occur. Our results 
corroborate the experimental observation, but in addi- 
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tion, we find that , at high temperatures, the ripples first 
come to rest before they are completely wiped out by the 
increasing non-linear contributions. 

The rest of the paper is organized as follows. First, 
we state our simulation model, i.e. how the sputtering 
process and the different diffusion mechanisms arc imple- 
mented. Then we explain, how we study the movement of 
the ripples. In the main section, we show our simulation 
results. We finish with our conclusions and an outlook. 



II. EROSION AND SURFACE MIGRATION 

According to Sigmund's sputtering theory 25], the 
rate at which material is removed from a solid surface, 
through the impact of energetic particles, is proportional 
to the power deposited there by the random slowing down 
of particles. The average energy E(r') deposited at sur- 
face point r' = (x r , y', —z') is given by the Gaussian dis- 
tribution 



E(r') 



(277)3/2(7^ 



■ exp 



(z' + d) 2 x' + y' 



2a 2 



(1) 



where we have used the local Cartesian coordinate sys- 
tem of the ion with origin at the point of penetration and 
with the z axis coinciding with the ion beam direction; 
[z 1 + d) is the distance of the surface point, from final 
stopping point of ion, measured along the ion trajectory, 
V x' 2 + y' 2 is the distance perpendicular to it; a and fj, 
are the widths of the distribution parallel and perpen- 
dicular to the ion trajectory respectively; e is the total 
energy deposited, d is the average depth of energy depo- 
sition. Sigmund's formula is the basis for all theoretical 
treatments and analysis of experimental results so far. 



A. The Sputtering Process 

Following [2H , we simulate the sputtering process on a 
surface of size L 2 with periodic boundary conditions, by 
starting an ion at a random position in a plane parallel 
to the plane of the initially flat surface, and projecting it 
along a straight trajectory inclined at angle 9 to the nor- 
mal to the average surface configuration; at an azimuthal 
angle 4>. The ion penetrates the solid through a depth d 
and releases its energy, such that an atom at a position 
r = (x, y, h) is eroded (see Fig. 1 of |24[) with probability 
proportional to E(r). It should be noted that, consistent 
with the assumptions of the theoretical models |2l l20ll2ll |. 
this sputtering model assumes no evaporation, no rede- 
position of eroded material, no preferential sputtering of 
surface material at point of penetration, and surface is 
defined by a single valued, discrete time dependent height 
function h(x,y,t) (solid-on-solid model, SOS). The time 
t is measured in terms of the ion fluence; i.e, number of 
incident ions per two-dimensional lattice site (x,y). Wc 
used incidence angle 9 = 50°, azimuthal angle <f> = 22.0°, 



d = 6.0, a = 3.0, fj, = 1.5, as obtained by SRIM (2(j 
for 5keV Xe + ions on graphite and rescaling all lengths 
by a factor 2. This should give according to the linear 
theory of Bradley and Harper a value 9 C = 68° |l2| . We 
have chosen e to be (27r) 3 / 2 cr/i 2 , which leads to high sput- 
tering yields Y ~ 7.0 compared to experiments like jll| . 
where Y = 0.3, ...,0.5, i.e. increasing the efficiency of 
the simulation. According to the Bradley Harper theory, 
the ripple wavelength A scales like A ~ Y~ x l 2 so that 
we expect patterns with correspondingly smaller length 
scales in our simulations. This we have to remember 
when quantitatively interpreting the result. Anyway, the 
general phenomena observed in the simulation are not 
affected by this choice. 

Our model of the sputtering mechanism sets the time 
scale of the simulation, and allows comparison with ex- 
periments. Additionally, also moves of atoms mimicking 
surface diffusion are performed, described now. 



B. The Hamiltonian and Arrhenius Models of 
Surface Diffusion 

Surface migration is modelled as a thermally activated 
nearest neighbor hopping process, as in |23,|28|]. A Monte 
Carlo acceptance/rejection procedure is used for this pur- 
pose. One diffusion step refers to a complete sweep of 
the lattice. Two different solid-on-solid models of sur- 
face diffusion in molecular beam epitaxy are used; the 
second one of them sensitive to the repulsion of a diffus- 
ing particle from a down step, and preferential diffusion 
in the uphill direction: the so-called Schwoebcl effect. 

The first model j2]j is based on a thermodynamic in- 
terpretation of the diffusion process. For each step, a site 
i and one neighbour site j are randomly selected. The 
trial move is an atom hopping from i to j, i.e. hi =hi — l 
and hj — hj + 1. We calculate the surface energy before 
and after the hop, through the energy of an unrestricted 
SOS model 



\hi - hjl 2 



(2) 



<*,]> 



J is a coupling constant through which the nearest neigh- 
bor sites interact, hi is the height variable at site i, and 
the summation is over the nearest neighbors on the 2- 
dimensional substrate. 

The hop is allowed with the probability 



pi-/ = V 



1 



exp h^ 



(3) 



where AE^f is the energy difference between the ini- 
tial and final states of the move. T is the substrate 
temperature, and ks is the Boltzmann's constant. Al- 
though no exact mapping is possible we can estimate that 
a temperature k^T/J ~ 0.2 in this model corresponds 
roughly to the temperature used in the second model be- 
low. The estimate is based on a comparison of the pure 
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diffusion mechanism without sputtering such that they 
fead to comparable values of the roughness. Note that 
this temperature is below the roughening transition of 
this model [27I • This model does not prevent atoms from 
moving down over steps edges, hence no Schwoebel effect 
is present. 

The second model is also based on a MC procedure and 
uses a formula known from kinetic MC mechanisms. For 
each step, again a site i and a nearest neighbor site j are 
chosen at random but now a hopping move is performed 
with a probability proportional to the hopping rate of an 
Arrhcnius form 



k(E, T) = k exp - 



— ) 

k B T 



(4) 



E = EsB+nnE^N + Es is an energy barrier to hopping, 
consisting of a Schwoebel barrier term Esb, a substrate 
term Eg = 0.75eV and a nearest neighbor bonding of 
magnitude nnE^N = un0.18eV; where nn is the number 
of in-plane nearest neighbors of the diffusing atom. Esb 
is equal to some constant (0.1 5eV in this case. Note 
that we perform also runs for Esb = 0, to compare 
with the thermodynamic model, see below), if the num- 
bers of next-nearest neighbors in the plane beneath the 
hopping atom before (nnrib) and after (nnn a ) the hop, 
obey 4 = nnni, > nnn a ; and zero otherwise. Our tem- 
perature is measured in units of eVfc^ 1 in this model, 
where T ~ 0.02eV^fc i3 1 corresponds to room temperature. 
ka = 2ksT/h is the vibrational frequency of a surface 
adatom, i.e. a hopping attempt rate, h being Planck's 
constant. The hopping attempt rate is very high, with 
a corresponding low hopping probability resulting from 
Eq. 01 slowing down the simulation. Thus we incorpo- 
rate the factor exp(— Es/kBT) into the rescaled attempt 
rate such that the hopping rate reads 



k(E,T) = fciexp 



AE 



(5) 



where k\ = fc exp 



Es 



is a much lower hopping at- 



tempt rate, AE = nnE^^ + EsB- This physical attempt 
rate, in comparison with the ion current density used in 
experiments, determines the ratio between the number of 
sputtering steps and the number of surface diffusion steps 
made in the simulation. In the next section, we state the 
values we used for our simulations. A discussion of pa- 
rameter optimization and a rcscaling of the temperature 
with the parameters is given in [29j Note finally that for 
atoms on top of planes, which are far from down edges, 
AE = 0, i.e. each hop is accepted, independently of the 
temperature. 



III. RIPPLE KINEMATICS 

In experiments, we typically have N = 
lxl0 15 atoms/cm 2 on the surface. Since the typi- 
cal experimental ion current density is of the order 



F = 7.5xl0 14 ions/cm 2 s [23|, this implies a flux of 
<t> = F/N ~ 0.75 ion/atom s. From the values given 
above, we get hopping attempt rates k\ of around 
200 1/s for room temperature, hence 200 sweeps of 
the diffusion mechanism correspond to 0.75 ions per 
surface atom. Thus, we initiate a diffusion step every 
$L 2 /fci = 0.0037L 2 erosion steps; L is the system size. 

Initially, for times less than about 1.4 ions/lattice site, 
the surface is rough [24| and then the formation of rip- 
ples starts. In this paper we focus on the motion and 
time development of these ripples. In Fig. QJthe time de- 
velopment of a sample surface topology is shown for the 
first diffusion model. Initially ripples are formed. They 
propagate slowly and, due to the increasing influence of 
non-linear effects (note the scales at the right), disappear 
at longer times. The long-time behavior, where the rip- 
ples have disappeared, has already been studied in Ref. 

m 




FIG. 1: Surface profiles at a substrate temperature of O^Jkg 1 
and at different times. Starting from top-bottom, left-right, 
t=0.5, 1.5, 4.0, 9.0, 14.0 and 20.0 ions/atom. Ion beam direc- 
tion, indicated by the bar, is perpendicular to ripple orienta- 
tion. The scales show the surface height measured from the 
lowest height. 

In order to monitor the ripple propagation on the com- 
puter, we assign the crest points of the ripples to clusters, 
and then monitor the motion of these clusters. A clus- 
ter of crest points is defined as the set of surface points 
with height h(x,y,t) > h c and nearest neighbor distance 
I < l c , where h c and l c are cut-off surface height and 
distance between neighboring cluster points respectively. 
We have chosen our cut-off height to be a function of the 
average height of the configuration (h), and the height 
difference hd between the maxima and minima of the sur- 
face; i.e, h c = (h)-\-phd, where p is some fixed percentage. 
In this way clusters with about the same proportionate 
sizes can be followed from the beginning of ripple for- 
mation until complete disappearance of the ripples. Fur- 
thermore, we have used l c = 2. Different, unconnected 
ripples should, in general, generate different clusters. We 
also require that the number N of elements in a cluster 
be large enough to allow for statistical analysis, here we 



4 



have chosen N > 10 elements. 



IV. RESULTS AND DISCUSSION 




FIG. 2: (Color online) Surface profile for time 3 ions/atom 
(thermodynamic diffusion model, T = 0.2J/fc_g, L = 128). 
The second profile contains the clusters formed from the first 
profile, as described in the text. 



The propagation of the ripples is studied by calculating 
the time rate of change of the position of the centre of 
mass of a cluster 



(6) 



where the summation is over all the elements of the clus- 
ter. We have assumed a homogeneous system composed 
of unit mass particles, such that the center of mass of a 
cluster is Xcm = Ei x i- The ripple wavelength is 
given by A = 2%/ri, r\ being the average expectation value 
of the Gaussian fitted to the peak of the structure factor 
S(k) = |/i(k)| 2 , where h(k) is the fourier transform of 
the height topography h(r,t), given by 



The results are obtained, as already mentioned, for 
square lattices of size 512x512, with periodic boundary 
conditions, and as an average over fifty different realiza- 
tions. 

For the case of the Arrhenius diffusion mechanism, (in- 
cluding the Schwoebel barriers) one can choose a temper- 
ature corresponding to the physical temperature present 
in the experimental system. A naive guess is to use room 
temperature fc^T = 0.02eV, at which the experiments 
usually are carried through. The resulting structures are 
shown in Fig. [3| for intermediate as well after long sput- 




FIG. 3: Sample surface topology for a small system (L=128), 
for the Arrhenius MC diffusion mechanism at surface temper- 
ature equal to room temperature, after t = 100 ions/atom. No 
clear ripples can be observed. Similar results were observed 
for almost all time, except the very early ones. 



fc ( k ) = ^7aE[ h ( r '*) 



(h)]e 



kr 



(7) 



d! is the substrate dimension, i.e. here d! = 2. Fig. [21 
shows two profiles of the surface for system size 128x128 
at time t — 3 ions/atom; in the second profile, we print 
the clusters on top of their corresponding ripples. As 
seen in the figure of the clusters, application of periodic 
boundary conditions neccessitates the need to first unfold 
toroidal clusters before calculating the position of their 
center of mass. As time increases, local surface slopes 
V/i increase, and since the non-linear effects depend on 
the square of \7h they will dominate by scaling down 
surface relaxation mechanisms |l| . These non-linear ef- 
fects are responsible for the disapearance of ripples (Fig. 
^) at long times, and for the transition of the surface 
topography from a periodic ripplepattern to a rough to- 
pography with self-affine scaling 0, Ha . We thus expect 
fluctuations in the position of the centre of mass due to 
disappearing ripples; the fluctuations are averaged out 
by using systems of size 512x512 with a large number of 
clusters such that the ripple velocity at any time is an 
average of the velocities of all the ripples at this time. 



tering times. We cannot observe clean ripples. The rea- 
son is that this kind of diffusion mechanism is too slow 
at room temperature to effectively counteract the strong 
roughening due to our model of sputtering, which pos- 
sesses a particularly high sputtering yield. Hops are al- 
most always prevented if an atom has in-plane neighbors, 
so the mechanism is not very effective on a rough surface. 
Since the surface relaxation is essential for the formation 
of ripples pfl |. it needs locally higher than room tem- 
peratures to produce clean ripples in our model. This 
happens indeed in experiments, since most of the kinetic 
energy, carried by the incoming ion, is converted into 
lattice vibrations, hence the surface is locally strongly 
heated. Here, we do not know the spatio-temporal dis- 
tribution of the local temperature. Either one would have 
to perform MD simulations, or include heat conduction 
in the model, both making the treatment of large systems 
over long time scales infeasible. Instead, we are choosing 
a higher but constant effective temperature T, which is 
a good first approximation. 

Now, we want to estimate this effective temperature. 
The most basic approach is to describe the energy carried 
by the ions as a constant inflow of energy at the surface, 
fix the temperature far away from the surface to room 
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temperature and solve the stationary heat-conduction 
equation to calculate the temperature at the surface [3(| • 
The resulting temperature depends strongly on the ion 
energy, the ion current density, and the thermal conduc- 
tivity of the material. For experimentally reported pa- 
rameters, temperature rises up to 1500 K (0.155eVfc^ 1 ) 
are found [3(j. This shows that high efficient temper- 
atures, even in the stationary state, may be achieved. 
However, in the experiments of Habenicht et al. |2^| 
only small average ion current densities have been used, 
which result in a temperature rise at the surface of only 
few K. 

This does not mean that one can use a temperature 
close to room temperature as effective temperature. The 
reason is that right after impact, the surface is strongly 
heated close to the melting temperature and the quickly 
cooled again, i.e. a thermal spike occurs |3l| . Further- 
more, the surface is sputtered using a focused ion beam 
(of diameter 30nm), which is moved relatively slowly 
over the surface and which exhibits a large spot cur- 
rent of 15/iA/cm 2 . Hence, under the ion beam, for sev- 
eral short time intervals, surface diffusion is greatly en- 
hanced. Marks has calculated [13 the spatio-temporal 
development of the temperature after ion impact by solv- 
ing the dynamic heat-conduction equation, resulting in a 
temperature profile T(r,t) as function of time t and dis- 
tance r from the point of impact. The initial distribution 
T(r, 0) is given by a step function with T(r, 0) being the 
melting temperature of the material for r < r$ and being 
the room temperature elsewhere. The initial radius ro is 
determined such that the thermal energy inside this semi 
sphere equals to the energy carried by the ion. Marks 
found that the surface is heated strongly right after the 
impact and is cooled down to temperatures close to room 
temperature within few ps. Qualitatively and quantita- 




FIG. 4: Surface profiles at a substrate temperature of 
O.leV/fcs with the second diffusion model. Starting from 
top-bottom, left-right, t=0.5, 1.5, 4.0, 9.0, 14.0 and 20.0 
ions/atom. In both cases, depicted here and in Fig. 1, ripples 
propagate along a direction opposite to that of the ion beam. 

tively similar profiles have been observed in MD simu- 
lations |3^| as well. We apply his equation, using the 



parameters for ion energy and ion current density in the 
spot as given above, to determine an effective tempera- 
ture (with ro = 15. 6 A in our case). The basic idea is 
that in an time interval At, the number of hops governd 
by the temperature T(0, t) at the impact point should be 
the same as under the effective temperature T: 

(8) 

We have neglected here the temperature dependence of 
fco- When including it, we found that the resulting effec- 
tive temperature changes only slightly. We have chosen 
At, as the average time between two ions arriving in a 
circle with area ttvq under the ion beam spot, resulting in 
At = 1.4 x 10 5 ps. For the energy barrier, we have chosen 
AE = E S b + 3E NN + Eg, which corresponds to atoms 
along edges of islands/steps. Using these parameters, we 
found an effective surface temperature of T — 1200K, 
i.e. considerably higher than room temperature. In this 
calculation it is assumed that only the energy carried by 
the ions hitting the "target area" tttq contribute to the 
heating of the surface inside the area. If one takes into 
account that also the ions hitting the neighborhood of 
the target area contribute to the heating inside the area, 
even higher effective temperatures can expected. 

The exact effective temperature depends on many pa- 
rameters as ion energy, ion current density, heat conduc- 
tion, surface roughness etc. We are here interested only 
in universal effects, not in modeling a specific experi- 
mental setup. For this reason, we use the above result 
only as a guideline and study several temperatures of 
this order of magnitude and additionally above. Hence, 
for the further analysis of ripple movement, we consider 
high effective temperatures for the Arrhenius MC model, 
such that the surface diffusion is indeed able to act as 
an effective smoothing mechanism (see Fig. 0J. At such 
higher temperatures we observe some universal features 
for both diffusion mechanisms, as presented now. Fig- 
urc[S]is the plot of the ripple wavelength (circle symbols) 
versus time measured in units of the number of ions per 
atom; its inset is a plot of the projection of the ripple 
velocity along the ion beam direction, versus time, both 
at the estimated effective temperature of fc^T = O.leV 
corresponding to the experimental conditions from Ref. 

m 

A plot of wavelength versus time in Fig. reveals that 
for short times A ~ i - 32 , which is in-between the results 
A ~ t 5 of Habenicht et al. and A ~ t 26 of Frost et 
al. |35|| . But we observed a power-law behavior only 
in the initial stages of ripple formation, the wavelength 
becoming constant in time at the later stage. 

The velocity shows a power-law behavior over a larger 
time interval, resulting in v ~ £~ - 7 as obtained from 
inset of Fig. [5] This is in excellent agreement with the 
experimental result v ~ t~ - 75 of Habenicht et al. [23| . 
A difference is that for smaller times a constant velocity 
was observed in the experiments, while we do not see any 
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t (ions/atom) 

FIG. 5: Ripple wavelength, A, measured in lattice units, as 
a fuction of time, t. The inset shows the time dependence of 
the ripple propagation velocity, v (measured in lattice units 
per ion per atom). Both results are for the kinetic diffusion 
mechanism, at a substrate temperature of feT = O.leV. 
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FIG. 7: Same plot as in Fig. |S|but for the Arrhenius diffusion 
mechanism, for a substrate temperature of 0.2 eV/kB- In both 
figures, the line with circle symbols represent the wavelength 
while the line with square symbols represent the velocity. 



clean ripples for smaller times than the power-law regime. 
Anyway, combining both scaling results gives v ~ A -219 , 
in good agreement with the exponent —2 of continuum 
theory Q. 




(ions/atom) 



FIG. 6: Ripple wavelength, A, measured in lattice units, as 
a fuction of time, t. The inset shows the time dependence of 
the ripple propagation velocity, v (measured in lattice units 
per ion per atom). Both results are for the thermodynamic 
diffusion mechanism, at a substrate temperature of 0.2 Jkg 1 . 



Now we turn to higher effective subtrate temperatures, 
corresponding e.g. to higher ion currents and/or materi- 
als with lower heat-conductivity. Figures and [7| are 
plots of the ripple wavelength (circle symbols) as a func- 
tion of time, at respective temperatures ksT = 0.2 J and 
kgT = 0.2eV; using the first and second models of sur- 
face diffusion respectively. In both models, the ripples 



disappear after a while, i.e. the ripple wavelength di- 
verges. Considering the lifetime of the ripples from first 
appearance to annihilation, the wavelength increases ex- 
ponentially with time as A ~ exp(pi), p = 0.029 (Fig|Hl), 
in the first model, while it increases with time according 
to the inverse law A(i) ~ 1/ (ci — c?£) with c\ = 0.083 and 
c2 = 0.0036 (Fig. 0) in the second model. To investigate 
the origin of the difference, we performed also simula- 
tions with the Arrhenius model, but with the Schwoebel 
term set to zero. In this case the result was very sim- 
ilar to result in Fig. HJ] of the thermodynamic model 
(which has no Schwoebel term here), and we obtained a 
behavior A ~ exp(0.036£). On the other hand, when we 
set the energy in the Schwoebel term to twice its value, 
Esb = 0.3eV, the result is very similar to E$b = 0.15eV. 
This shows that the Schwoebel barrier plays an important 
role in the pattern formation process. Note that these 
fits are purely heuristic. We are not aware of any theory 
of the time dependence of ripple wavelength and veloc- 
ity, only a calculation of the dispersion relation v(A) has 
been performed within linear theory Q. Furthermore, 
there exists an analytic study of the temporal develop- 
ment of step bunches during epitaxial growth [23, 13 • 

The insets of Figs. ED and [7| are plots of the ripple ve- 
locity (line with square symbols) as a function of time. 
Irrespective of which surface diffusion mechanism is em- 
ployed, the velocity is at first almost independent of time, 
then it disperses after a transition time t r . This initial 
plateau is similar to the plateau observed in the experi- 
ments, but the drop in velocity is very abrupt, no clear 
power law is visible then. Moreover, the ripples finally 
come to rest before completely disappearing, as seen in 
the smaller inset of Fig. HD We find, however, that at 
the lower temperature in the kinetic model, the ripples 
do not stop moving until their disappearance. Figure [S] 
shows the dependence of the ripple velocity on the wave- 
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FIG. 8: Ripple velocity as a function of ripple wavelength, 
for the thermodyanmic and, in the inset, for the Arrhenius 
surface diffusion mechanism. 



length for ksT — 0.2eV resp 0.2J, their order of magni- 
tude relationship is about the same as in the experiment. 
We see in Fig. [5] that the trend in velocity variation is 
the same at high temperatures but the magnitude in- 
creases with temperature, as one would expect from the 
temperature dependence of the surface diffusion. But we 
only observed a power-law scaling at temperatures be- 
low fceT ps 0.18eV. This indicates that the presence of 
power-law scaling of ripple wavelength and velocity, and 
the corresponding exponents, depend on the time scale 
of observation (Fig. |SJ, as well as on the effective tem- 
perature. 

It seems that the increase in magnitude of the velocity, 
when measured at same time (t < t r ) but different tem- 
peratures, does not continue indefinitely in our model. 
In the upper graph of Fig. there is very little differ- 
ence in the magnitudes of the velocity at temperatures 
2.0 and 5-OJfc^ 1 ; even though the temperature difference 
is very high. This saturation behaviour is also displayed 
in the ripple wavelength at the same higher tempera- 
tures, as seen in the lower graph of Fig. |5J In principle, 
one can still fit an exponential law to the data, except 
that the decay constant p in the exponential becomes 
very small (It has the respective values of 0.029, 0.018, 
0.0031, and 0.003 from the lowest to the highest temper- 
ature.). So for very high effective temperatures we could 
equally well fit a power-law. Hence, there may be some 
"critical substrate temperature" , above which the wave- 
length remains nearly constant in time; and the velocity, 
after some time t r , drops instantaneously to zero. Nev- 
ertheless, the temperature where such a "transition" will 
take place, is probably unphysically high (see below), so 
that the material used in the experiment would start to 
evaporate before reaching this point. But other materi- 
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FIG. 9: Temperature dependence of the ripple velocity (up- 
per graph), and ripple wavelength (lower graph); for the ther- 
modynamic diffusion mechanism. Temperature is in units of 



als, in combination with high ion currents, might quite 
display such a behavior. So far, we are only aware of one 
set of experiments p3j . Hence, it would be very inter- 
esting to see, whether some temperature dependence of 
the dynamical features, including the disappearance of 
the coarsening, can be seen in experiments at higher ef- 
fective temperatures corresponding to high ion currents 
and/or higher lab temperatures. 

Our results for the second diffusion model also indicate 
that in X(t) « l/(ci - C2^), C2 approaches zero with in- 
creasing temperature. Here, where we can measure the 
temperature in real units, it is clear that the "transition" 
to almost non-coarsening ripples, takes place at unrealis- 
tically high temperatures 2 — 5eV/ItB, where the material 
starts to evaporate. Moreover, we notice in Fig. [^jthat 
the transition time from linear regime to onset of non- 
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linearities decreases with increasing temperature. 

To summarize, ripple propagation depends on the ef- 
fective substrate temperature as well as diffusion mecha- 
nism. At around so-far experimentally realized tempera- 
tures, ripples propagate, from first appearance, with de- 
creasing velocity until disappearance without full cessa- 
tion of motion. At high effective temperatures, however, 
immediately after ripple formation, the ripples move with 
constant velocity for some time, after which they begin 
to decelerate (insets of Figs. El Q and after some time, 
depending on the diffusion model, the ripples stop mov- 
ing but keep disappearing gradually. At the same time 
the ripple structure is gradually being washed out, and 
in the final stage the ripples are completely wiped out. 
The ripple wavelength is always increasing in time at 
high temperatures, while at low effective temperatures it 
initially increases with time, and later becomes constant. 

V. CONCLUSION AND OUTLOOK 

We have studied the propagation of ripples by means 
of a discrete (2 + l)-dimensional model of the sputtering 
process, combined with one of two different solid-on-solid 
models of surface diffusion: a Arrhenius MC mechanism 
with Ehrlich-Schwoebel barriers and a thermodynamic 
mechanism without a Schwoebel term. We have obtained 
the formation and propagation of the ripples with both 
diffusion mechanisms used in turn. Furthermore, we have 
obtained the same trend in the behavior of ripple veloc- 
ity and wavelength as observed experimentally and pre- 
dicted theoretically, but, in addition to the experimental 
results, we find a drastic change in the ripple propaga- 
tion at temperatures well above the so-far experimentally 
realized effective temperature; for instance we found de- 
viations from power-law into exponential or inverse-law 
behavior, and in addition, the ripples first stop moving 
before vanishing completely. We find that, at very high 
effective temperatures, the behavior of the ripple velocity 
is charaterized by two regions, separated at the transition 
time. In the first region it is constant and in the second 
region it decreases rapidly to zero. Between the two re- 
gions a power-law dependence can be observed for some 
small time interval. Whereas, around so-far experimen- 
tally realized temperature, the velocity-time relationship 
obeys a power-law. Furthermore, at high effective tem- 
peratures, the wavelength increases exponentially with 
time in the thermodynamic diffusion model (and in the 
Arrhenius diffusion model without Schwoebel term) and 
obeys an inverse law for the Arrhenius model including 
the Schwoebel barrier. In addition, we find further strong 
dependencies on the effective substrate temperature; as 
the temperature increases the magnitude of the veloc- 



ity also increases. The transition time between constant 
and decreasing velocity is also found to decrease with in- 
creasing temperature. Our results indicate an approach 
towards a saturation behaviour of velocity or wavelength 
with increasing effective substrate temperature, where 
the wavelength is expected to become time independent. 
However, this may happen at an unphysically high tem- 
perature. Anyway, an experimental study of the depen- 
dance of the dynamical features of ripple formation and 
effacement on the physical conditions seems very promis- 
ing. 

One open problem of our model, at high incidence an- 
gles (e.g 6 > 75°), is that it uses the Sigmund formula for 
modeling the sputtering process. In a recent simulation 
[3^ | using a binary collision approximation, we observed 
that close to the penetration point of the ion, much less 
atoms were sputtered than predicted by the Sigmund for- 
mula . in fact the distribution shows a minimum there. 
When incorporating this effect in the Bradley-Harper lin- 
ear theory [2(j of sputtering, one e.g. observes |3(| that 
the sputter yield, i.e. the number of removed atoms per 
ion, exhibits a minimum for grazing incidence, like in ex- 
periments, in contrast to the orginial linear theory |20| . 
Hence, it may be promising to apply a different formula 
describing the sputtering, which takes into account this 
effect. 

Furthermore, the role of the interplay between the sur- 
face diffusion process and the sputtering process is still 
not fully understood. So far, we know that including 
a pure T = relaxation in our sputtering model does 
not [24j lead to a disappearance of ripples for long times. 
Next, we know from this study, that one approach includ- 
ing calibrated Schwoebel barriers does not yield ripples at 
room temperature for sputtering yield Y ~ 7. There are 
several different models [13, HE HE HE ISE E3, E3, E3| for 
surface diffusion, which could be combined in a construc- 
tion kit manner. Here, an extensive study over different 
combinations of parameters is necessary. 

Finally, it would be of interest to include crystal 
anisotropy into the surface diffusion; This may give re- 
sults in agreement with experimental studies of metallic 
substrates, which may be useful in understanding the 
anomalies of such surfaces. 
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